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Abstract. The dynamical equations describing the evolution of a self-gravitating 
fluid of cold dark matter (CDM) can be written in the form of a Schrodinger equation 
coupled to a Poisson equation describing Newtonian gravity. It has recently been shown 
that, in the quasi-linear regime, the Schrodinger equation can be reduced to the exactly 
solvable free-particle Schrodinger equation. The free-particle Schrodinger equation 
forms the basis of a new approximation scheme -the free-particle approximation - that 
is capable of evolving cosmological density perturbations into the quasi-linear regime. 
The free-particle approximation is essentially an alternative to the adhesion model 
in which the artificial viscosity term in Burgers' equation is replaced by a non-linear 
term known as the quantum pressure. Simple one-dimensional tests of the free-particle 
method have yielded encouraging results. In this paper we comprehensively test the 
free-particle approximation in a more cosmologically relevant scenario by appealing 
to an A'^-body simulation. We compare our results with those obtained from two 
established methods: the linearized fluid approach and the Zeldovich approximation. 
We find that the free-particle approximation comprehensively out-performs both of 
these approximation schemes in all tests carried out and thus provides another useful 
analytical tool for studying structure formation on cosmological scales. 
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1. Introduction 

Over the last two decades a vigorous interplay between theory and observation has 
led to the establishment of a standard cosmological model that not only describes the 
composition and evolution of the universe as a whole, but has also scored great successes 
in accounting for the large-scale structures observed in the universe; for a recent review 
see, for example, [1]. In this standard model, large-scale structure is thought to be 
the result of the gravitational amplification of small density perturbations present in 
the early universe; a mechanism known as gravitational instability. The gravitational 
instability process is relatively easy to understand at a qualitative level because it mostly 
relies only on the action of gravity, with more complex hydrodynamical and radiative 
effects playing a role only at late times when structure is already largely developed. 
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Moreover, when the density fluctuations are very much smaller than the average density, 
the growth of fluctuations can be handled quite accurately using first-order (linear) 
perturbative techniques. On the other hand, galaxies and galaxy clusters represent 
density enhancements well in excess of the mean density so we need to go beyond linear 
theory if we are to study the formation of such objects in detail. 

The non-linear regime of gravitational instability is analytically intractable, except 
in a few cases where some special symmetry applies. For the most part, therefore, 
cosmologists have been forced to resort to numerical methods based on A^-body 
computations in order to study the late stages of the evolution of density fluctuations. 
The use of computer simulations has revolutionized cosmology, especially in enabling the 
simulation of the large-scale structure expected to be found in galaxy redshift surveys 
and the consequent testing of models using observations. In some cases, however, such 
as when a very large volume and/or very high resolution is needed, one has to make 
more detailed predictions than can be obtained by using computational techniques. In 
such cases numerous analytical approximations have been suggested that provide robust 
results in a necessarily limited range of circumstances (e.g. see [21 [3] and references 
therein). These approximate techniques also help to promote a genuine understanding 
of the principal processes involved in the formation of cosmological structure. 

Analytical attempts to understand large-scale structure formation are commonly 
based on the assumption that cold dark matter (CDM) can be treated as a self- 
gravitating pressureless fluid. Two important approximation methods can be derived 
from the equations of motion for this fluid by applying linear Eulerian and Lagrangian 
perturbation theory, respectively. First-order Eulerian perturbation theory - the so- 
called linearized fluid approach (e.g. [4J) - has been the backbone of structure formation 
theory for many years as a result of its simplicity and robustness at early times/on 
large scales. However, for Gaussian initial conditions, it generically leads to unphysical 
regions where the matter density is negative if extrapolated to late times/small 
scales. Although popular, the linearized fluid approach is considerably less powerful 
than its Lagrangian counterpart - the Zeldovich approximation [5]. The Zeldovich 
approximation considers (linear) perturbations in the trajectories of individual fluid 
elements (particles), rather than in macroscopic fluid quantities such as the density 
field. This technique is remarkably successful in describing the pattern of structure 
formed from realistic (random) initial conditions (e.g. [6j), as long as care is taken to 
ensure that particle trajectories do not intersect. Once trajectories cross (a phenomenon 
known as shell- crossing), the velocity field becomes multi- valued and the density field 
develops a singularity (known as a caustic). However, since the Zeldovich approximation 
is kinematical, particles do not respond to the strong gravitational forces acting in the 
vicinity of the caustic. Instead they simply continue to move along their original inertial 
trajectories and any structure formed is rapidly 'smeared out' in an unrealistic manner. 
An alternative to the Zeldovich approximation which avoids this problem is the adhesion 
model [7], an ad hoc extension of the Zeldovich approximation in which particles are 
assumed to 'stick' to each other when shell-crossing occurs. This sticking is achieved by 
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including an artificial viscosity term in the usual hydrodynamical equations of motion, 
thus transforming the momentum equation into Burgers' equation. The viscous term 
becomes large when particle trajectories intersect, cancelling the component of the 
particle velocities perpendicular to the caustic. This prevents particles from simply 
sailing through the caustic and, consequently, the density singularities predicted by the 
Zeldovich approximation are regularized and stable structures are formed. Comparisons 
with A^-body simulations have shown that the adhesion model is capable of accurately 
reproducing the skeleton of the 'cosmic web' of filaments and voids that comprises 
the morphology of the large-scale galaxy distribution [U [U [TOl [TTl [12] . Although the 
original motivation for the adhesion approximation was purely phenomenological, it 
has recently been shown that Burgers' equation can be naturally derived from the 
coarse-grained equations of motion, provided certain simplifying assumptions are made 
[I31[T1[I51[I61[I7]. 

A novel alternative approach to the study of large-scale structure formation was 
suggested by Widrow and Kaiser [1^. They proposed a wave-mechanical description 
of self-gravitating matter in which coUisionless CDM is modelled by a complex scalar 
field whose dynamics are governed by coupled Schrodinger and Poisson equations (see 
[T9] for a relativistic extension of the original theory) . As pointed out by Coles [20] , the 
(approximately) log-normal form of the distribution function of density perturbations 
can be naturally explained within this wave-mechanical formalism. By appealing to 
simple one-dimensional examples of gravitational collapse in a static universe. Coles and 
Spencer [21] showed that the wave-mechanical approach offers a competitive alternative 
to the Zeldovich approximation, provided the effective Planck constant u in the theory is 
chosen carefully. They also demonstrated that, unlike the Zeldovich approximation, the 
wave-mechanical approach leads to a density field that remains non-singular when shell- 
crossing occurs. Building on the work of Coles and Spencer, Short and Coles [22] recently 
showed that, in an expanding universe, the Schrodinger equation can be conveniently 
reduced to the free-particle Schrodinger equation in the linear and quasi-linear regimes 
of gravitational instability. They advocated using the free-particle Schrodinger equation 
as the basis of a new approximation scheme - the free-particle approximation. The free- 
particle approximation is essentially an alternative to the adhesion model in which the 
viscosity term is replaced by a non-linear term known as the quantum pressure; it is this 
term that prevents the formation of caustics. A convenient property of the free-particle 
Schrodinger equation is that, like Burgers' equation, it possesses an analytic solution. 
Short and Coles [22] used the free-particle method to follow the gravitational evolution 
of a plane-symmetric sinusoidal density fluctuation until shell-crossing. By considering 
such a simple example they were able to carefully elucidate the effect of the quantum 
pressure term. It was found that, if not properly controlled, the quantum pressure term 
significantly impedes the growth of density perturbations. However, in the case where 
the effect of this unusual term was minimized, the free-particle approximation provided 
an excellent match to the Zeldovich approximation (which is exact in one-dimension 
up until particle trajectories intersect) and was capable of generating large (non- linear) 
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over-densities. Our objective in this paper is to provide the first comprehensive test of 
this promising new method in a cosmologically relevant scenario by appeahng to a full 
A^-body simulation. We also assess the performance of the free-particle approximation 
relative to the established linearized fluid approach and the Zeldovich approximation. 

The layout of the paper is as follows: In section [2] we present the central equations 
of the free-particle approximation and comment on the link between the free-particle, 
adhesion and Zeldovich approximations. The essentials of the practical implementation 
of the free-particle method are discussed in section [3l In section H] we describe the results 
of our comparison of the free-particle approximation with an A^-body simulation. We 
draw our conclusions in section [5l 



2. The free-particle approximation 

In the free-particle approximation, coUisionless CDM is represented by a complex scalar 
field ip = D) obeying the free-particle Schrodinger equation 

where the amplitude and phase of the wavefunction are related via 

^V^[argW] + -^(|^r-l)=0, (2) 

and the effective Planck constant i/ is a free (real) parameter with dimensions of L^. For 
a full derivation of these equations and a discussion of the validity of the free-particle 
approximation, we refer the reader to [22]. The spatial coordinates x = x(D) are 
comoving coordinates and we use the linear growth factor D = D(t) a,s a, 'time' variable, 
rather than cosmological proper time t. The linear growth factor is the growing mode 
solution of the differential equation 

b + 2Hb - ATiGpb^cD = 0, (3) 

where a dot denotes a derivative with respect to t and, at some initial time ti, 
D[ = D{ti) = 1. The Hubble parameter H = H{t) is defined by if = a/a, where the 
scale factor a = a(t) is normalized to unity at the present epoch, and pb.c = Ph,c(t) is the 
density of CDM in a homogeneous Friedmann- Robertson- Walker background cosmology. 
The wavefunction is assumed to be of the Madelung form [23] 

^ = (1 + 5)1/2 expf^Y (4) 



where the density contrast 5 = (5(x, D) is defined hj 6 = p/ph,c — 1 and describes 
fluctuations in the CDM density field p = p(x, D) about the homogeneous background 
value pb,c- The gradient of the velocity potential = 0(x, D) gives the comoving velocity 
field u = u(x, D), defined by u = dx/dD. Inserting the Madelung transformation (jlj) 
into and ([2]) yields 
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-Vx-[(l + 5)Vx0]=O 



(6) 



(5) 



and 



5 = 



(7) 



respectively. The term V = ^(x, D) is known as the quantum pressure term and is 
given by 



The hnk between the free-particle approximation and the more conventional Zeldovich 
and adhesion approximations becomes apparent upon inspecting the Bernoulli-like 
equation If we were to set V = then reduces to the so-called Zeldovich- 
Bernoulli equation [24j which, along with the relation ([7D, arises when the Zeldovich 
approximation is expressed in Eulerian, rather than Lagrangian, space. On the other 
hand, replacing P by a term of the form /iVx0 leads to Burgers' equation [25] for the case 
of an irrotational velocity field; this is the defining equation of the adhesion model (e.g. 
[7]). The term /xVx'/' prevents the formation of multi-stream regions and regularizes 
the density singularities predicted by the Zeldovich approximation. It is controlled 
by the (real) parameter /i which formally plays the role of a viscosity coefficient in 
the adhesion approximation. The growth of density perturbations is suppressed on 
scales < /i^/^ by the viscous term and thus, for large values of /i, the adhesion model 
leads to a much smoother distribution of matter than observed in A^-body simulations 
(e.g. [9j). In order to maximize the dynamic range of the adhesion approximation, 
the (inviscid) limit ;U — > is commonly employed (e.g. [3, El [HI EE]). In this limit, 
the adhesion model reduces exactly to the Zeldovich approximation away from regions 
where particle trajectories intersect. In modern approximation methods (such as the 
Euler- Jeans-Newton model [I3] and the small-size expansion [151 IE]) where Burgers' 
equation is consistently derived (rather than assumed), the constant viscosity coefficient 
/i is replaced by a density-dependent gravitational multi-stream coefficient (see [E] and 
references therein for a detailed discussion). 

The quantum pressure term V cannot be written in a form proportional to V^^ 
and so the free-particle and adhesion approximations are not equivalent. However, 
the quantum pressure term does play a qualitatively similar role to the viscous term 
in the sense that it is also a regularizing term, preventing the formation of multi- 
stream regions and caustics. This was demonstrated by Coles and Spencer [21] who 
used the Schrodinger equation to follow the evolution of a one- dimensional sinusoidal 
density fluctuation beyond shell-crossing. The parameter v controls the quantum 
pressure term in the same way that /x controls the viscous term and, like /z, v can 



V 



2 (1 + 5)1/2 



(8) 
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also be viewed as an approximation to a general gravitational multi-stream coefficient. 
Short and Coles [22] showed that, in the semi-classical limit u 0, the free-particle 
approximation reduces exactly to the Zeldovich approximation prior to shell-crossing. 
The quantum pressure term then only becomes important in regions where particle 
trajectories intersect. However, for finite u, the quantum pressure term has an effect 
before shell-crossing occurs and acts to inhibit the gravitational collapse of density 
perturbations. The larger the value of u, the greater the suppression effect; see [22] 
for a discussion. The free-particle approximation is then no longer identical to the 
Zeldovich approximation, although they may become similar at a certain distance from 
collapsing regions, depending on the actual value of z/. Also, once shell-crossing occurs, 
the quantum pressure term has an effect outside of multi-stream regions in the finite u 
case (cf. figure [T] of f2Tj). 

3. Testing the free-particle approximation 

In modern cosmology the process of large-scale structure formation is almost exclusively 
studied using numerical A^-body methods since these techniques allow the gravitational 
evolution of density perturbations to be followed far into the non-linear regime. This is 
something that cannot be done with any analytical approximation. As a result, A^-body 
simulations provide an ideal tool for testing the applicability of analytical approximation 
schemes. In this paper we test the free-particle approximation, the linearized fluid 
approach and the Zeldovich approximation in Eulerian space (which we henceforth refer 
to as the Zeldovich-Bernoulli approximation) against an A^-body simulation so that we 
may assess the performance of the new free-particle approximation relative to existing 
methods. 

3.1. The N-body simulation 

The adaptive P^M code HYDRA [27] was used to perform an A^-body simulation. 
The simulation consisted of A^ = 128^ CDM particles contained within a cubic box 
of comoving side length L = 200h~^ Mpc equipped with periodic boundary conditions. 
Here h is the present value of the Hubble parameter Hq in units of 100 km s^^ Mpc~^. 
A spatially-flat ACDM cosmological model was assumed with parameters Qc,o = 0.24, 
^b,o = 0, ^^A,o = 0.76, h = 0.73 where flc,o, ^b,o and f^A^o are the present values of the 
CDM, baryon and dark energy density parameters, respectively. The simulation began 
at an initial scale factor Oi = 0.02 and the particles were initially positioned so as to 
form a random realization of a CDM density field with a power spectrum P = P{k) of 
the form: 

P = T^Pi, (9) 

where Pi = Pi{k) is the primordial power spectrum of density fluctuations generated 
by inflation and T = T{k) is the transfer function. We assume that the primordial 
power spectrum is scale- invariant: Pi oc k, as predicted by a general class of inflationary 
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models (see, for example, [2B] for a review). The following empirical form for the CDM 
transfer function was adopted [29] : 

T = ^Wl^^"^^ [1 + 3.89g + (16.1g)^ + (5.46g)^ + (6.71g)^] "^^^ , (10) 

where q = q{k) is given by g = kO^/"^ / {VLcoh'^ Mpc~^) and 6 = k/1.68. Here n = 
i7r^o/^7,o is the ratio of the present relativistic particle and photon density parameters 
respectively; in a universe with three relativistic neutrino flavours and photons k, = 1.68. 
Finally, the power spectrum was normalized so that the present value of the rms density 
fluctuation in spheres of radius 8h~^ Mpc (crg^o) was 0.74. The Newtonian gravitational 
evolution of the CDM distribution was then followed from ai = 0.02 to the present 
Oo = 1, with the comoving positions of all N particles stored at 20 different output 
values of the scale factor. 



3.2. The testing process 

To test the Eulerian free-particle, linearized fluid and Zeldovich-Bernoulli approxima- 
tions against the Lagrangian A^-body simulation described above, we divide the cubic 
simulation volume (of comoving side length L = 200h~^ Mpc) into a uniform cubic 
mesh with A^^g = 128^ grid points. There is then one particle per grid cell (on av- 
erage) and the comoving separation between neighbouring grid points in the mesh is 
A = L/Ng^^ = 1.5625/1-1 Mpc, which is roughly the size of a typical galaxy group. At 
each output value of the scale factor, we use the comoving positions of the particles 
in the simulation to generate the CDM density field on the cubic mesh by applying 
triangular-shaped cloud interpolation (e.g. [30]). The three approximation methods we 
are comparing all use the linear growth factor D as a time variable, rather than the 
scale factor a. We therefore convert the output values of the scale factor used in the 
simulation to values of the linear growth factor by noting that, in a spatially- flat ACDM 
cosmology, the solution of ([3]) is 

where the constant of proportionality is chosen to ensure Di = 1, Ba is the incomplete 
beta function and a = a{a) is defined by 



^ = 7\ — — ZE- (12) 



We now describe how the CDM density field can be calculated at any particular value 
of the linear growth factor for each of the approximation schemes we are considering. 
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The free-particle approximation. The first step is to determine the initial velocity 
potential (pi on the cubic mesh. This is done by numerically solving V^^i = Si in 
Fourier space, where Si is the initial CDM density field of the A^-body simulation. The 
initial wavefunction ipi is generated on the mesh by inserting Si and 0i into 

^i = (l + 50^/'expf^Y (13) 



along with a suitable finite value of u. We discuss how the optimal value of u can be found 
in section I4T1 The discrete Fourier transform ipi = ipi(k) of the initial wavefunction is 
then calculated and the exact solution of the free-particle Schrodinger equation ([T]): 



%[) = tpi exp 



(14) 



is used to determine the Fourier transform 'ip = '?/'(k, D) of the wavefunction at any 
D > 1 of interest. Here k is a co moving wavevector and k = |k|. The wavefunction ip 
in real space follows by taking the inverse discrete Fourier transform of ip and the CDM 
density field is determined from the amplitude of the wavefunction via S = {ipl"^ — 1. 
As we shall see in section HI knowledge of the velocity potential (p is also useful. We 
calculate (p on the cubic mesh by using numerical Fourier techniques to solve ([Tj). 



The linearized fluid approach. The CDM density field S at any D > 1 is determined 
from the initial A^-body density field Si by simply applying the linear growth law S = DSi. 



The Zeldovich- Bernoulli approximation. The initial velocity potential (pi is constructed 
on the cubic mesh as before and the Zeldovich-Bernoulli equation 

|^-i|V«^|-0 (15) 

is integrated forwards in time from = 1 to Dq using the the simple numerical scheme 
described in [21]. At each time step in the integration we know (p and so the CDM 
density field S can be numerically calculated on the mesh via ([7]). 



4. Results and discussion 



In order to optimize the free-particle approximation, it is essential to select the 
parameter v carefully. We briefiy address this issue before proceeding to test the free- 
particle method against the linearized fiuid and Zeldovich-Bernoulli approximations. 
First, we introduce a dimensionless parameter F = z//A^ with A defined as before. The 
parameter F emerges upon rewriting the free-particle Schrodinger equation ([1]) in terms 
of the dimensionless comoving coordinate x = x/A. Note that F oc i^; hereafter we will 
analyse the behaviour of the free-particle approximation in terms of F rather than v. 



Wave mechanics and the adhesion approximation 



9 



4-1- Optimizing the free-particle approximation 

The numerical implementation of the free-particle approximation described in section 
13.21 requires a finite value of F, rather than F — 0. Recall that, in this case, the 
quantum pressure term impedes the gravitational collapse of density perturbations even 
before multi-stream regions are formed. However, since P oc F^, this undesirable effect 
can be minimized by choosing the smallest possible value of F. In order to find this 
optimal value of F, we first note that the phase of the initial wavefunction ( fT3|) is 
arg (-i/^i) = — 0i/z/ oc 1/F. Samphng the initial phase field at the Nyquist sampling 
rate requires that the change in phase between two adjacent grid points in the cubic 
mesh must be less than or equal to vr radians. As F is decreased towards zero, the 
phase of the initial wavefunction varies increasingly rapidly and this condition will 
eventually be violated. Phase-aliasing effects then cause the power spectrum of the 
initial wavefunction ip^ to become very noisy and our method breaks down. The optimal 
value of F, Fc, is then defined as the smallest value of F for which there is no phase- 
aliasing; we ensure there is no aliasing by using a simple numerical algorithm to check 
that the phase difference between each point in the mesh and its nearest neighbours 
is not greater than tt radians. The value of Fc is mesh- dependent and, in the work 
presented here, F^ = 0.12. To illustrate that F = Fc is indeed the optimal case, we 
compare the relative sizes of the quantum pressure V and convective C = — | Vx0|^ /2 
terms in the Bernoulli equation ([5]) for several different choices of F. We do this by 
introducing a ratio x = x{(^)y defined by x = For each value of F we also 

investigate how the quantum pressure term affects the formation of large-scale structure 
within the free-particle approximation. 

Figure [1] shows how the ratio x evolves with the scale factor. For the moment 
we focus on the three solid lines, corresponding to F = 0.12, F = 0.6 and F = 20, 
respectively. Observe that the ratio x increases monotonically with the scale factor 
in all cases, implying that, on average, the quantum pressure grows relative to the 
convective term as time proceeds. We can also see that, at any particular value of 
the scale factor, the value of x decreases as F is decreased. Therefore, on average, 
the quantum pressure term becomes less significant (relative to the convective term) 
as the value of F is made smaller. At each value of the scale factor, we find that the 
smallest value of x is indeed obtained for F = Fc = 0.12. For reasons that will become 
clear in section 14.21 we also examine the behaviour of x when the free-particle density 
field is smoothed with a Gaussian filter prior to calculating the quantum pressure and 
convective terms (note that we apply the smoothing first to ensure numerical gradients 
are well behaved). The dashed and dotted lines in figure [1] show how x evolves for two 
smoothing scales: = 4:h~^ Mpc and rgm = 8h^^ Mpc, respectively. In both cases the 
value of X is very small at all values of the scale factor and, as the smoothing length 
is increased, the quantum pressure term diminishes relative to the convective term (on 
average). 

In figure [2] we show slices through the final (i.e. at a = oq = 1) A^-body and free- 
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Figure 1. The ratio x — plotted as a function of the scale factor a for 

three different values of the parameter F. The dashed and dotted lines show how x 
evolves when the free-particle (with F = 0.12) density field has been smoothed with 
Gaussian filters of radii = Ah~^ Mpc and rgm — 8h~^ Mpc, respectively. 

particle density fields for the cases F = 0.12, F = 0.6 and F = 20. The global morphology 
of the free-particle density field agrees well with that of the iV-body density field for 
F = 0.12, with the peaks of the two density fields located at the same positions. However, 
the over-densities formed in the free-particle approximation are not as large as in the 
A^-body simulation and the free-particle density field is considerably smoother than the 
A^-body density field. This is not surprising since the free-particle approximation is 
essentially a quasi-linear method and thus cannot account for the full non-linear growth 
of density perturbations. As the value of F is increased to 0.6 the pattern of large-scale 
structure appears 'smeared out' in the free-particle approximation. This is because the 
quantum pressure term is now larger on average (see figure [1]) and inhibits the collapse 
of density perturbations over a greater distance; see also [22]. Similar behaviour is also 
seen in the adhesion approximation as the viscosity parameter fi is increased [9]. In the 
case where F = 20, the quantum pressure term completely dominates the convective 
term on average (figure [I]) and there is no growth of density fluctuations whatsoever. 
Instead, the initial density field simply oscillates rapidly in a seemingly random fashion 
and remains close to homogeneous. This oscillatory behaviour was also observed in the 
simple test of the free-particle approximation performed by Short and Coles [22] . 
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Figure 2. Slices through the A^-body density field S^h and the free-particle density 
fields 6ip for three different choices of the dimensionless parameter F. The value of 
the scale factor is a = oq = 1 and the slices are arbitrarily taken at a z-coordinate of 
lOOh-^ Mpc. 

4-2. Comparison of approximation schemes 

We have established that our free-particle method will be optimized for F = Fc = 0.12. 
Accordingly, we set F = 0.12 in the free-particle approximation for the remainder of 
this work. We now use the method described in section 13.21 to test the free-particle 
approximation against an A^-body simulation and compare our results with those of the 
linearized fluid approach and the Zeldovich-BernouUi approximation. 

Figure [3] shows point-by-point comparisons of the density fields obtained from the 
three approximation schemes with the A^-body density field. If the density field 5 
obtained from any one of our approximation methods is identical to the A^-body density 
field (5nb at all points in the cubic mesh, then there will be no scatter about the diagonal 
line overlaid on each plot. We quantify any scatter about the diagonal by using the 
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Figure 3. Point-by-point comparisons of the density fields obtained from the three 
approximation schemes and the A^-body density field (5nb- The first, second and third 
columns correspond to the free-particle (with F = 0.12), linearized fluid and Zeldovich- 
BernouUi approximations, respectively. Comparisons are shown for three different 
values of the scale factor a and all of the density fields have been smoothed with a 
Gaussian filter of radius = 4/i~^ Mpc. The corresponding value of the correlation 
coefficient rg between the density fields is shown in the top- left corner of each plot. 



correlation coefficient 

(5^)1/2(52)1/2' l^"^^ 

a value of = 1 corresponds to the case where there is no scatter. All of the density 
fields in figure [3] have been smoothed with a Gaussian filter of radius = 4:h~^ Mpc 
to smooth over highly non-linear regions where the linear /quasi-linear approximation 
methods we are testing cannot be expected to give reliable results. The main result of 
figure [3] is that the free-particle approximation performs significantly better than both 
the linearized fiuid and Zeldovich-Bernoulli approximations when tested against an A^- 
body simulation. It is clear from the point-by-point comparisons that, at any given 
value of the scale factor, the free-particle density field displays the best correlation with 
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the A^-body density field. Observe that the extrapolation of first-order Eulerian and 
Lagrangian perturbation theory to late times leads to unphysical regions in the density 
field where S < —1, i.e. where p < 0. This problem is neatly side-stepped in the free- 
particle approximation since, by definition, S = {ipl"^ — 1 and lipl"^ > 0. This explains 
why the free-particle density field provides the best match to the A^-body density field in 
low density regions. It is evident from figure [3] that, even after applying smoothing, all 
of the approximation methods systematically underestimate the density in high-density 
regions. However, it is interesting to note that the free-particle approximation performs 
considerably better than the other approximation methods in high density regions too. 
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Figure 4. Point-by-point comparisons of the density fields obtained from the three 
approximation schemes and the A^-body density field 5nb- The layout of the plots is 
the same as in figure HJ the only difference is that all of the density fields have been 
smoothed with a Gaussian filter of radius = 8/i~^ Mpc. 



Figure H] is essentially the same as figure [3] except that all of the density fields 
have been smoothed with a Gaussian filter of radius = 8h'^^ Mpc. By using a 
larger smoothing scale we restrict all of the density fields to the quasi-linear regime 
5 ~ 1. Consequently, the performance of all the approximation methods is improved. 
However, the free-particle approximation still clearly out-performs the linearized fluid 
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and Zeldovich-BernouUi approximations. The correlation between the smoothed free- 
particle and A^-body density fields is very good at all values of the scale factor, with only 
a slight discrepancy between the two density fields in high-density regions as a — >■ 1. 
Although the smoothing is now much heavier, the linearized fiuid approach and the 
Zeldovich-Bernoulli approximation still generate negative matter densities at late times. 

The correlation coefficient provides a simple way of quantifying the scatter about 
the diagonal in the point-by-point comparisons of figure [3] and figure |H In figure [5] we 
explicitly plot rs as a function of the scale factor for each approximation scheme. Results 
are shown for the two smoothing scales: r^^ = 4h~^ Mpc and = 8h^^ Mpc. It is 
immediately clear that, for both smoothing lengths, the free-particle approximation 
consistently yields a higher value of rs than the other two methods. It is interesting to 
note that, when the free-particle density field is smoothed with a Gaussian filter of radius 

= 4:h~^ Mpc, the value of rs obtained is always larger than that calculated from 
the linearized fiuid and Zeldovich-Bernoulli density fields smoothed on a scale of 8h~^ 
Mpc. In other words, even when the linearized fiuid and Zeldovich-Bernoulli density 
fields are more heavily smoothed, the correlation with the A^-body density field is still 
not as good. 
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Figure 5. The correlation coefficient rs for each of the three approximation schemes 
plotted as a function of the scale factor a. Results are shown for the two smoothing 
lengths — 4/i^^ Mpc and rsm — 8h^^ Mpc. 



A simple means of examining the statistical properties of the density fields in the 
different approximation schemes is afforded by the one-point probability distribution 
function (PDF) Ps = PsiS) of the density fields. The PDF of the initial A^-body density 
field is Gaussian by construction. However, non-linear gravitational evolution causes the 
A^-body PDF to become non-Gaussian at late times. The distribution becomes skewed 
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since positive density fluctuations (over-densities) can grow indefinitely whereas negative 
density fluctuations (under-densities) cannot exceed 6 = —1. Figure [6] shows the PDFs 
of the various density fields at the final time a = oq = 1 for the smoothing lengths 
'"sm = 4:h~^ Mpc and = 8h~^ Mpc. The corresponding PDFs of the iV-body density 
field are shown for reference. We can see that, for both smoothing scales, the PDF of the 
free-particle density field provides the best match to the A^-body PDF for all values of S. 
The fact that the free-particle approximation guarantees 6 > 1 causes the free-particle 
PDF to be skewed in the same way as the A^-body PDF. It is clear from both plots in 
figure [6] that the PDF of the linearized fluid density field is Gaussian. This must be the 
case for Gaussian initial conditions since the density contrast at a given point simply 
grows proportional to the linear growth factor to first-order in Eulerian perturbation 
theory. The PDFs of the Zeldovich-Bernoulli density field also seem to be close to 
Gaussian. Consequently, the linearized fluid and Zeldovich-Bernoulli approximations 
both assign a non-zero probability to the existence of regions with 5 < — 1 at late times. 
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Figure 6. The one-point PDFs Pf, of the density fields obtained from the free-particle 
(with r = 0.12), linearized fluid and Zeldovich-Bernoulli approximations. The A''- 
body PDFs are also shown for comparative purposes. The value of the scale factor 
is a = Oq = 1 and the top and bottom plots correspond to the smoothing lengths 
Tsm = 4/i~^ Mpc and fsm = 8/i~^ Mpc, respectively. 

In our discussion so far we have only compared the performance of the free-particle, 
linearized fluid and Zeldovich-Bernoulli approximations in real space. However, it is 
informative to examine the behaviour of the different approximation schemes in Fourier 
space as well. For each method we investigate how the Fourier components b = 5(k, a) 
of the appropriate density field evolve relative to the corresponding Fourier components 
^nh = '^nb(k, a) of the A^-body density field. To do this we use the statistic e = e{k, a), 
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Figure 7. The statistic e plotted as a function of tlie comoving wavenumber k for 
each of the three approximation schemes. The value of the scale factor is a = oq = 1. 
The dashed vertical line corresponds to the wavenumber fc„i = 0.3927/i Mpc~^ which 
indicates the boundary between Fourier modes in the linear and non-linear regimes. 

defined by 

^-'^ . (17) 

where, for eacfi (comoving) wavenumber fc, tfie summations are over all wavevectors k 
with a magnitude in the interval (k — 2tt/ L,k\. The quantity e provides a measure of 
the differences between the amplitudes and phases of the Fourier components 5 and 5nb- 
If the amplitudes and phases of the Fourier components of the two fields are identical 
then e = 0. On the contrary, if the phases of the two fields are completely uncorrelated 
then e = 1 on average. A convenient property of the statistic e is that it is independent 
of any smoothing of the density fields and so we are free to use the un-smoothed density 
fields to calculate e. Figure [7] shows e as a function of the wavenumber k at the final time 
a = Oq = 1 for the different approximation schemes we are considering. The vertical 
dashed line shown in the plot denotes the wavenumber fcni = 0.3927/i Mpc~^ which 
corresponds to a wavelength of Aa;„i = 16/i~^ Mpc. From observations, the present rms 
density fiuctuation in spheres of diameter Afc^j is of order unity and so the scale /cni can 
be taken as an indication of the boundary between modes in the linear and non-linear 
regimes of gravitational evolution. Fourier modes with a wavenumber k < k^\ are still 
in the linear regime and evolve independently of each other. The phase information 
of the initial density field is preserved in this case. On the other hand, modes with a 
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wavenumber k > kj^i experience phase shifts due to couphng between different modes 
in the non-hnear regime (e.g. [311 E^)- apparent from figure [7] that, for all of 

the approximation schemes, the value of e becomes larger as k is increased. This is not 
surprising since the various methods we are comparing cannot accurately follow the non- 
linear evolution of small-scale modes. The main point of figure[7]is that, for wavenumbers 
in the range < < 3A;nb the value of e is consistently smaller in the free-particle 
approximation than in the linearized fluid and Zeldovich-Bernoulli approximations. In 
particular, the free-particle approximation performs much better than the other two 
methods for < A; < k^i. Notice that the Zeldovich-Bernoulli approximation gives a 
slightly lower value of e than the free-particle approximation for k > 3fcni- However, 
the wavenumber k = 3k^\ corresponds to a wavelength Xsk^i = 5.3333h~^ Mpc and we 
find that the rms density fluctuation in spheres of diameter X^k^^ is approximately 2 for 
the final A^-body density field. Thus we cannot realistically expect the free-particle and 
Zeldovich-Bernoulli approximations to yield reliable results on scales much smaller than 
this anyway. 

5. Conclusion 

In this paper we have investigated a new approach to the study of cosmological structure 
formation known as the free-particle approximation. The free-particle approximation is 
similar in spirit to the adhesion approximation in the sense that it is an extension of 
the Zeldovich approximation which includes a regularizing term - the quantum pressure 
term - to prevent the formation of density singularities when shell-crossing occurs. The 
quantum pressure is controlled by a free parameter u (or F in the notation of section 
H]) which is akin to the viscosity coefficient fi in the adhesion model. In the semi- 
classical limit u —>■ the quantum pressure only has an effect in multi-stream regions. 
However, in this work we have employed a numerical implementation of the free-particle 
approximation that uses a finite value of u instead. In this case the quantum pressure 
term acts to suppress the collapse of density perturbations away from regions where 
particle trajectories cross; the degree of suppression increases with u. 

The free-particle method has been rigorously tested by appealing to a full 
cosmological A^-body simulation. We have shown that the performance of the free- 
particle approximation, relative to the A^-body simulation, is optimized when u is set 
to the smallest possible value supported by the discrete mesh we use. In this case 
the effect of the quantum pressure is minimal. For larger values of u > Uc the quantum 
pressure term becomes more significant relative to the convective term in ([5]) and the 
details of the large-scale structure distribution appear 'washed out' in the free-particle 
approximation. In particular, for z/ 3> z/c, the quantum pressure term is dominant and 
initial density fluctuations undergo oscillation rather than growth. The free-particle 
approximation (with u = Uc) consistently yields better results than the established 
linearized fluid and Zeldovich-Bernoulli approximations in all of the tests discussed 
in section HI At all times the density field obtained from the free-particle method is 
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considerably better correlated with the A^-body density field than the density fields of 
the other two approximation schemes. This is particularly true in low density regions 
which is a consequence of the fact that the free-particle density field is guaranteed to 
be non-negative. This feature of the free-particle approximation also explains why the 
one-point PDF of the free-particle density field provides the closest match to the A^- 
body PDF of all the approximation methods. We have also seen that the free-particle 
approximation out-performs the linearized fluid approach and the Zeldovich-BernouUi 
approximation in Fourier space. From a practical point of view, another advantage of 
the free-particle approximation is that it is very quick to implement. This is because the 
free-particle Schrodinger equation has an exact solution, allowing us to directly calculate 
the density and velocity fields at any given time without having to numerically integrate 
equations of motion. 

To conclude, the free-particle approximation provides a quick and effective way of 
evolving cosmological density perturbations into the quasi-linear regime in any CDM- 
dominated cosmology. We have compared the free-particle method with two traditional 
approaches and the results have proved favourable. In light of these comments, 
we believe that the free-particle approximation provides another useful addition to 
the repertoire of analytical techniques available for the study of large-scale structure 
formation. 
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